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A new approximate solution to the quantum-classical Liouville equation is derived 
starting from the formal solution of this equation in forward-backward form. The 
time evolution of a mixed quantum-classical system described by this equation is 
obtained in a coherent state basis using the mapping representation, which expresses 
N quantum degrees of freedom in a 2A'"-dimensional phase space. The solution yields 
a simple non-Hamiltonian dynamics in which a set of coherent state coordinates 
evolve in forward and backward trajectories while the bath coordinates evolve under 
the influence of the mean potential that depends on these forward and backward tra- 
jectories. It is shown that the solution satisfies the differential form of the quantum- 
classical Liouville equation exactly. Relations to other mixed quantum-classical and 
semi-classical schemes are discussed. 



I. INTRODUCTION 



Nonadiabatic processes are at the core of many physical phenomena, including population 
transfer among electronic system states, quantum coherent evolution of a system interacting 
with environmental degrees of freedom, electron and proton transfer reactions in condensed 
phase and biological systems, among others. In investigating such phenomena one often 
focuses on certain quantum degrees of freedom whose dynamics is of primary interest. These 
may be the electronic degrees of freedom of a chromophore excited by radiation to prepare 
the initial state of the system, the exciton states of a light harvesting system, or even the 
electron or proton degrees of freedom involved in the transfer of these particles. In such cases 
we are led to consider how these quantum degrees of freedom interact with the environment 
in which they reside. Interactions with the environment can lead to the breakdown of the 
Born-Oppenheimer approximation and one must consider nonadiabatic dynamics in such 
open quantum systems. 

A number of different approaches have been developed to describe nonadiabatic dynam- 
ics. These include mean-field and a variety of surface-hopping scheme^^, methods based 
on semi-classical evaluations of path integral formulations of quantum mechanicJ^^ and 
descriptions based on the quantum-classical Liouville equation!^. An important ingredi- 
ent in any approach dealing with nonadiabatic dynamics is the manner in which quantum 
coherence and decoherence are taken into account in the dynamics. The description of nona- 
diabtic dynamics necessarily entails dealing with coherence that is generated and destroyed 
as the system evolves while interacting with its environment. Many of the various nonadi- 
abatic approaches that have been constructed deal with the issue of decoherence in various 
waypHSSl. 

Another characteristic of nonadiabatic schemes is the manner in which the environment 
is modeled. At the simplest level, the environment may be treated as a stochastic bath, 
which leads to reduced descriptions that do not explicitly include the environmental degrees 
of freedom in the evolution. Their effect only appears in certain parameters and terms that 
characterize the coupling to the environment. Schemes of this type include various quan- 
tum master equations=^^, the Lindblad equatiorP^ and the Redfield and Bloch equationd^HMl 
Other methods explicitly account for the environmental degrees of freedom. It is challenging 
to treat large and complex systems fully quantum mechanically, although there are devel- 
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opments along these lineJ^^^. Some methods, for example, some path integral methods, 
begin with a full quantum treatment and then make semi-classical approximations to obtain 
tractable solution Often the environment in which the quantum dynamics of inter- 
est occurs can be described by classical dynamics to a high degree of accuracy and this has 
spawned a number of mixed quantum-classical descriptions of nonadiabatic dynamics. Many 
surface-hopping schemes fall in this category as do some approximations to semi-classical 
path integral formulations and mean-field methodj ^^ * ^^ * ^^ l Here we focus on descriptions 
based on the quantum-classical Liouville equation (QCLE). 

The QCLE employs a partial Wigner representation of the environmental (bath) degrees 
of freedom and may be derived from full quantum dynamics by truncating the quantum evo- 
lution operator to first order in a small parameter related to the ratio of the characteristic 
masses of quantum and bath degrees of freedomP^. It may also be derived from partially 
linearized path integral formulationd^QEi]^ indicating the close connection between these dif- 
ferent starting points. This equation has been shown to provide an accurate description of 
nonadiabtic dynamics in many applications and to account for quantum decoherencd^. A 
number of different methods, whose structure depends on the basis chosen to represent the 
quantum degrees of freedom, have been devised for its simulationH^HlSl. Simulation methods 
that utilize an adiabatic basis can be cast into the form of surface-hopping dynamics, but in 
a way that includes coherent evolution segments that account for creation and destruction 
of coherence in a proper manner. More recently, as in some semi-classical approached, the 
mapping basi^ was used to describe the quantum degrees of freedom in the QCLE in a 

continuous classical-like manner, leading to a trajectory description in the full system phase 
spac^SHU 

In this paper we also utilize the mapping representation but instead of dealing directly 
with the solution of the QCLE using a Liouville propagator, we start with its solution 
in terms of forward-backward quantum-classical propagators constructed some time agd^. 
With this starting point and the introduction of a coherent state basiJ^we are able to obtain 
a solution of the QCLE that involves forward-backward trajectories of the coherent state 
variables, coupled to the evolution of the bath phase space variables. Formally, both forward 
and backward trajectories are propagated forward in time. The two sets of trajectories are 
distinguished and named by their association with the forward and backward quantum- 
classical propagators, respectively. This formulation leads to a simple set of non-Hamiltonian 
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equations that describe the nonadiabatic dynamics of the system. 

The outhne of the paper is as follows: In Sec. [TT] we sketch the important features of 
the QCLE, its representation in the mapping basis and formal solution in forward-backward 
form needed for our calculation. The forward-backward trajectory solution is constructed 
in Sec. |III[ which contains the most important results of the paper. A discussion of the 



results is presented in Sec. |IV[ while the Appendices give additional technical details of the 
calculation. 



II. QUANTUM-CLASSICAL LIOUVILLE EQUATION 

We consider a quantum subsystem coupled to a bath. We assume that the dynamics 
of such a system is described by the quantum-classical Liouville equatiorP^'^^SSESlMI 
a quantum operator By\/{X), which depends on the classical phase space variables X = 
{R, P) = {Ri, R2, ...,Rn^, Pi, P2, -PatJ of the bath, this evolution equation takes the form, 

j^Bw{X,t) = iCBw{X,t), (1) 

where the quantum-classical Liouville operator is 

tC- = '-[Hw, ■] - li{Hw, ■} - {-, Hw}). (2) 

Here the subscript W refers to a partial Wigner transform over the bath degrees of freedom 
(DOF), H\y{X) is the partial Wigner transform of the total Hamiltonian of the system, 
[■, ■] is the commutator and {■, ■} is the Poisson bracket in the phase space of the classical 
variables X. The total Hamiltonian may be written as the sum of bath, subsystem and 
coupling terms, 

Hw{X) = H,{X) + hs + V,{R) = H,{X) + h{R), (3) 

where Hb{X) = P^/2M + Vb{R) is the bath Hamihonian with Vb{R) the bath potential 
energy, hg = /2m + Vs is the subsystem Hamiltonian with p and Vs the subsystem momen- 
tum and potential energy operators, and Vc{R) is the couphng potential energy operator. 
The masses of the subsystem and bath particles are m and M, respectively. The evolution 
equation for the density matrix pwiX,t) is analogous to Eq. ([T| with a change in sign of 
the evolution operator. 
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A. Formal solution 



The QCLE may also be written in a form that is analogous to the quantum Liouville 
equationlSS 

j^BwiX,t) = ^(^nABw-BwnA), (4) 
where operators "Ha and 'Ha are given by 

nA = Hwil + —j, nA=il + ^jHw, (5) 

<— — > <— — >■ 

with A the negative of the Poisson bracket operator, A =Vp ■ Vi?, — V_r ■ Vp- 
The formal solution of the QCLE can be expressed in either of two forms as 

Bw{X,t) = e'^'Bw{X) (6) 
= S {e'^^''''Bw{X)e~'^^''''^ 

The first equality follows from the formal solution of Eq. ([T]) while the second equality 
follows from Eq. Q. The S in this latter form simply specifies the order in which products 
of the left and right operators act in order to be identical with the first from involving 
the QCL operator. In particular, a general term S (^(j-LAy Bw{'Ha)''^ in the expansion of 
the exponential operators is composed of ^'^j^ separate terms each with a prefactor of 
Each of these separate terms corresponds to a specific order in which the "Ha and "Ha 
operators act on Bw- This formal solution will be used in the calculations presented below. 



B. Mapping representation 

We will be concerned with the representation of the QCLE in the quantum subsystem 
basis and its equivalent representation in the mapping basis. The subsystem basis, {|A); A = 
1, . . . , A^}, is defined by the eigenvalue problem hs\X) = €x\X), and a matrix element of an 
operator BwiX) is given by B^' (X) = {X\BwiX)\X'). 

The I A) eigenfunctions of an A^-state quantum subsystem can be replaced with eigenfunc- 
tions of A^ fictitious harmonic oscillator^^^ESl^ I^a), having occupation numbers which are 
limited to or 1: |A) — )■ Inix) = |0i, ■ ■ ■ , 1a, ■ ■ ■ Oat). Creation and annihilation operators on 
these states, d\ and dx, respectively, are defined as 



and satisfy the commutation relation [dx, dy] = (5a, a'- The actions of these operators on the 
single-excitation mapping states are d\ |0) = \mx) and d\ \mx) = |0), where |0) = |0i . . . On) 
is the ground state of the mapping basis. 

We may then define mapping versions of operators, Bm{X), given by 

BUX) = B^''{X)d{dy, (8) 

so that a matrix element of Bw in the subsystem basis is equal to the matrix element of 
the corresponding mapping operator in the mapping single-excitation basis: B^' (X) = 
{X\I3]y{X)\X') = {mx\Bm{X)\mx') . (The Einstein summation convention will used through- 
out although sometimes sums will be explicitly written if there is the possibility of confusion.) 
In particular, the mapping Hamiltonian operator is 

Hm = H,{X) + h^^\R)d\dx> = H,{X) + h^, (9) 



where we applied the mapping transformation only on the part of the Hamiltonian that 
involves the subsystem DOF in Eq. The pure bath term, Hb{X) in Eq. acts as an 
identity operator in the subsystem basis and is mapped onto the identity operator of the 
mapping space. 



The QCLE Q may now be written in terms of mapping operators as 

j^BUx, t)='-(^niB„,- Brn ) , (10) 
— ^ — >■ ^ 

where is given by = Hm{^ + ^A/2z), with an analogous definition for H™. One 
may verify that the mapping space matrix elements of this equation are identical to the 
subsystem matrix elements of Eq. Q. Consequently, the formal solution of this equation is 
similar to that in Eq. ^ and is given by 

BUX, t)=S (^e^«"X'*/^5„(x)e-^«"A */^^ . (n) 

This equation will form the starting point for the explicit solution of the QCLE in terms of 
forward-backward trajectories. 



III. FORWARD-BACKWARD TRAJECTORY SOLUTION 



The formal solution of the QCLE can be written in terms of a sequence of M short-time 
propagators acting on the initial value of the operator: 

Bw{X, t) = e^^^*ie*^^*2 . . . e'^^'^' Bw{X), (12) 

where Atj = tj — tj-i = t for all j with to = and Im = t. (When information about a 
specific time step is needed we use the Atj notation, otherwise the common value r will be 



used.) Consequently, in view of Eq. (11), the formal solution applies in each time segment 
so Bw(^X,t) may also be written as 

RUX, t)=S (e^^^^^^^A (e^^^^^^A . . . 



where there are M concatenated S (■ ■ ■) brackets. 



A. Representation in coherent states 

In order to proceed with the evaluation we must consider the computation of the forward 
and backward propagators in this expression. To order we have 

^irn^/n ^ ^H^AT/2^iHmT/h _^ o^t"^). (14) 

Also, to order we may write the first exponential operator as 

e^'"^^/^ = l + ^^^A + ..., (15) 
= 1 + ^H,{X)A + ^h'^^'alayA + . . . , 



1 + ^Ht,{X)A + ^ (^h^^'aya{ - Ti, A + 



where we have reversed the normal-ordered product of annihilation and creation opera- 
tors into an anti-normal order form using their commutation relation. The by-product of 
reversing the ordering of creation and annihilation operators is the emergence of a trace 
term in the last line of this equation. Since the trace term is independent of the quantum 



state, it may be combined with the bath potential, Vo(i?) = Vb{R) — Tis h{R), to give 
Ho{X) = P^/2M + Vo{R) so that we have the simpler form of Eq. (jisj). 



^ ^ ^ ^Ho{X)A + ^h^^'ayd[A + 0{t^). (16) 

In this form, the propagator can be expressed conveniently in coherent stated^. 
We define the coherent states 1^;) in the mapping space, 

ax \z) = zx \z) , {z\ a\ = zl {z\ , (17) 

where is a coherent state with degrees of freedom and the eigenvalue is zx = {qx + 
ipx)/V2h. The variables q = (gi, . . . , q^) and p = (pi, . . . ,p]\f) are mean coordinates and 
momenta of the harmonic oscillators in the state 1^;), respectively; i.e., we have {z\qx\z) = qx 
and {z\px \z) = px- 

The coherent states form an overcomplete basis; thus, we have to specify the inner product 
between any pair of coherent states and the resolution of identity^^. The inner product is 

= g-|(k-^T)-i3(^-2'*)_ (;^g) 

The norm of the inner product measures how far away the two coherent states l^;) and \z') 
are in the phase space of coherent state variables. The resolution of the identity is 



\z) {z\ , (19) 



where d^z = d{^{z))d{Q{z)) = dqdp/{2h)^. 

Given these properties of the coherent states, we may insert the resolution of the identity 
in the bath Hamiltonian terms and between the dy and operators in Eq. ( 16 ) to obtain 

e^^-^ = (1 + IHo{X)A) I ^ 1^) {z\ 

+ 1 1 ^h>^''^y\^) {z\d{A + 0{r^) 

= /5k) {l + lmX) + h'''zlzy)A 



'^''^■;2)e5^-(^'^)^(^| + 0(r2). (20) 



In this calculation we used Eq. (17) to eliminate the annihilation and creation operators in 



Eq. (20). Note that h'^'^' z\z\i = ^h'^'^' (qyqx + pxpy) since h'^'^' is symmetric. In the last line 



of Eq. ( 20 ) we defined the "classical" Hamiltonian 



H,iiX, z) = HoiX) + h^^ zlzy = HoiX) + ha{R. z) 

p2 



(21) 



where V,i{R,z) = Vo{R) + V,""^' {R)zlzy. 

The operator Hci{X, z)A acts on all bath phase space variables to its right. Therefore, it 
is convenient to introduce a notation that makes this action evident. More specifically, we 
let 



Hci{X, z}A = —— ■ — — ■ — =1 C {X,z) 



dP dR OR dP 



so that 



(Pz 



71 



N 



z) e 



iC{X,z)T/2 



Similarly we can define 



(22) 



(23) 



AH,i{X,z) 



d dH,i d OH, 



d 



dP OR dR dP 



-I L (X, z) 



(24) 



and 



d^Z 



N 



z) e 



iC{X,z)T/2 



z\ + 0{t^ 



(25) 



The other quantity that will enter in the evaluation of the time evolution is the action of 
the exponential operator Q^Hm{x)T/h ^ coherent state. In Appendices A and B we show 
that 



-iHm{X)T/n 



-iHt{X)T/h-ihm{R)r/n 



'iHt{X)T/h 



(26) 



with z{t) determined from the solution of the evolution equation. 



dz\ i dhr 



dt h dz\ 



(27) 
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B. Time evolution of an operator 



These results may now be used to compute the value of the matrix elements of an operator 
Bw{X,t) in the subsystem basis: By^'{X,t) = {mx\Bm{X,t)\my) . We have 



(41 e-*^™^ 14) e'^(^'^i)^) (41 mv), (28) 

We may now make use of the definition of the S operator to rewrite the actions of the 
right and left operators acting on the bath coordinates of an arbitrary operator A\y{X) in 
terms of a single effective operator C(.{X, z') that depends on the coherent state variables 
z and z' associated with the forward and backward propagators, respectively. In Appendix 
C we show that 

S (^e*2(^.-) i Aw (X)e*^(^'^') 5 ) (29) 
= e"^^^''^'^''^^Aw{X) = AwiXr). 

The explicit form of iCe{X, z^z') is 

■r (Y '\ P ^ dVe{X,z,z') d 



where Ve{X, z, z') = {Vd{R, z) + Vd{R, z'))/2. From Eqs. ([29]) and ([30]) we can see that the 
time evolution of the bath coordinates under the effective Liouville operator is given by the 
solutions of the equations 

dR^P_ dP^ dVeiX,z,z') 

dt M' dt dR ' ^ ^ 
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These results may be used in the expression for B^'{X,t) in Eq. (1281) to 



„ M 
iiii' i=l 



(mx \zi) {z[\mx') 



{zi\e 



Ati 



\Z2) 

■ ■ ■ \zm) 



(32) 



-iH„ 



14) ) (4|e 



-iH„ 



14) 



This expression can be evahiated by applying the operators from left to right. For example, 
the action of the first effective bath operator updates the bath phase space coordinates from 
X = XtQ to Xt-^ . Thus, 



„ M 



dPz' (P z'- 

:i]v^:i#("^A ki) (41 "^A') 



(zi|e*^™(^*i)^ \Z2) 
...S{^^'(X,J...(4|e-^^-(^*i)^|4) 



(33) 



The coherent state matrix elements can now be evaluated using Eq. (26) to give 

M 



nil' i=l 



{mx \zi) (41 "^A') 



s{^^'(x,j . . . 14) )e-^^(^*i)^*^/^ (41 4(tO) 

-Y. I ^^(^Akl)(4l"^A') 



(;^i(ti)k2)e^^=(^'-^-^^)^((;^2|... 

i?{^^'(x,j...|4))(4i4(tO)). 



(34) 



In writing the last equality we canceled the phase factors involving Hb{XtJ. 

At this point we can see how a description involving continuous trajectories may be con- 
structed. The classical bath propagator for the next time step from ti to ^2, e*^'"*^^'i'^^'^2)^^ 
involves the coherent state phase space variables Z2 and z'2 which may take any values 
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from the set of coherent state values. The coherent states involved in the matrix elements 
(zi(ti) 1^2) and (2:21 z[(ti)) are not orthogonal since the coherent states are overcomplete. 



However, in view of Eq. (18), we see that the overlap between two coherent states decays 
rapidly if their phase space coordinates differ significantly. Consequently we assume that 
(zi(ti) 1^2) ~ 7!'^6{z2 — zi{ti)) and (2:21 z[(ti)) ^ tt^6{z'2 — z[(ti)). Then performing the 
integrals over Z2 and Z2 we obtain 

B'^\X,t) = Y,j (35) 

2 = 3 



x[{z,{t,)\...B^^^\x,,)...\z[{t{)))). 



All coherent state and bath phase space variables have now been updated to time ti and 
process can now be repeated for all M time steps, starting with the application of the 
effective bath evolution operator for the time step At2- The result of this process is the 
simple expression 

B'^\X^) = Y^j ^^{mx ki) {z[\my) (36) 



xi^{z,{t)\m^)Bi:^' {Xt){m,,\z[{t)) 

The matrix elements between coherent states and the single-excitation mapping states may 
be evaluated explicitly to give 

(m,|z) = z,e-l^l'/^ (37) 

Writing this expression in terms of the x = {q,p) variables, and using the fact that "^^{ql + 
pI) is conserved under coherent state dynamics, we obtain 

B^'{X,t) = XI / dxdx'<p{x)(P{x') 

^^Mt)-^P,mq'At)+'Mt))^ (38) 

where 4>{x) = {27ih)~^ e~^''^'^'^~^^^^^'^^ is the normalized Gaussian distribution function and 
we have removed the subscript 1 from the dummy coherent state variables. The coupled 
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equations of motion governing this evolution are 

dqx ^ dHci{R,P,q,p) dpx ^ dHa{R,P,q,p) 
dt dpx ' dt dqx 

dq', dH4R,P,q',p') dp\ dH,i{R,P,q',p') 



dt dp'^ ' dt 

dR_P_ dP _ dHe{R,P,q,p,q',p') 
~dt ~ M' ~dt ~ M ' 



(39) 



where 



H,{R,P,q,p,q',p')= (40) 
^iH,i{R,P,q,p) + H,i{R,P,q',p')). 



Equation (38 ) and the associated evolution equations (39 ) are the results we set out to derive. 
They constitute a simple algorithm for obtaining a solution to the QCLE. Figure [T] presents 
a schematic picture that depicts the dynamics of coordinates prescribed by the evolution 



equations (39). As noted earlier, although both forward and backward trajectories are 
propagated forward in time, the two sets of trajectories arise from the forward and backward 
quantum-classical propagators, respectively. 

Earlier it was shown that the solution to the QCLE in the mapping basis can be given 
in terms of an ensemble of entangled trajectories. The solution in Eq. (39) is consistent 



with this interpretation in that the forward and backward trajectories of the coherent state 
variables are linked by the evolution of the bath variables and the evolution equations are 
in non-Hamiltonian form. A more detailed link between these two different approaches to 
the QCLE in the mapping basis is a topic that merits further study. 

C. Back to Differential Form 

In this section we show that the solution constructed above is indeed a solution of the 
QCLE. We do this by deriving the QCLE in the subsystem basis by constructing a finite- 
difference expression for the time evolution of By^' (X,t). We first write the matrix element 



for (A I Bw{X, t + r) |A') using Eq. (|36|), 

B^^\X,t + r) = J2l d^z{t)d^z'{t)<P{z)<P{z') (41) 
xz,{t)z':,{t)z;{t + r)z'^,{t + r) (e^^-i?{^^'(X, t)^ 
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where (f){z) — n 1^1^/^. We then expand to first order in r to obtain 

B^^'{X,t + T)^Yl [ d^z{t)d^z'{t)cf>{z)cf>{z') 



XZx"^ 



+r[z;it)^ + z',it)^)B^/iX,t) 

+Tz;{t)z'^,{t)zCeB^^'{X,t)\ + 0{t'). (42) 

The integrals over z{t) and z'{t) may be performed and, after rearranging terms and taking 
the hmit r ^ 0, the result is (some details are given in Appendix D), 

hm -B^ it) (43) 

^{X\'-[Hw,Bw] \X') 

~2 ^^^^ \^^w,Bw^ - ^Bw,Hw^ |A')^ , 

which is the QCLE. 

The QCLE in the subsystem basis is a first order differential equation with respect to time; 
therefore, it only describes how the matrix elements of B\y{X, t) at the beginning and the end 
of a time step are related. That our solution is found to satisfy the QCLE is consistent with 
the fact that all approximations used to derive the evolution in a single time step are exact 
to C(t^). However, in order to connect the trajectories of coherent state phase variables 
from adjacent time steps, we made the approximation, {zi{T) \ -Zi+i) ~ 7r^6{zi^i — Zi^T)). To 
understand the effects of this approximation, wc consider how our solution would be modified 
if the approximation were not made. One way to re-formulate the solution is to insert a 
set of single-excitation mapping states between every inner product of coherent states, i.e. 
{zi{T) \ Zi+i) — X]^. {zi{T) \ rUm) (m^.| -Zj+i). Once the mapping states are inserted, one loses 
the continuous trajectory picture in the coherent state phase space but one can formally 
integrate out the Zi and z[ variables in sequential (or chronological) order. This sequence of 
formal integrations is equivalent to evaluations of the matrix elements of 5^^/ (X, t) at every 
time step. Computationally, this is a very demanding task because one needs to sample, 
propagate and integrate out coherent state trajectories at every time step. However, this 
prescription (a continuous evolution of matrix elements) coincides exactly with the dynamics 
one would expect from the QCLE in the subsystem basis. 
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At this point, it is obvious that the coherent-state orthogonahty approximation re- 
places the continuous evolution of the matrix elements, B^' {X,t), with continuous tra- 
jectories, z{t) and z'{t). Instead of taking B!^ {X,t — r) as the starting point to compute 
B^'(X,t) at the next time step. The orthogonality approximation actually takes the oper- 
ator \z{t — r)) {z{0) \ Bm{X, 0) |2;'(0)) {z'(t — r) \ as the starting point and further propagates 
trajectories from the previous time step to obtain \z(t)) {z{0)\ Bm{X,0) |-2'(0)) {z'(t)\. Al- 
though the orthogonality approximation inevitably yields nonlocal errors, it does provide a 
computationally efficient way to simulate the dynamics. Other semi-classical approaches for 
solving the system-bath dynamics indicate that this is a sensible approximation to make. 
For instance, if we do not use the orthogonality approximation then we can write our solu- 
tion in the form of a standard coherent state path integral. Application of the stationary 
phase approximation will yield the same set of equations of motion for the coherent state 
phase variables. Similar coherent state dynamics was obtained in the context of a different 
semi-classical frameworlP^. 

Finally, we comment on the fact that the semiclassical analysis yields exact quantum 
mechanical solution for quadratic Hamiltonians. This is certainly true when the system is 
isolated from the bath. The same also holds true for our solution; if there are no bath terms 
then there is no need to make the orthogonality approximation. However, when a bath 
is present, the semi-classical analysis is equivalent to implicitly making the orthogonality 



approximation, which becomes exact in the limit /i — )■ in view of Eq. (18). The potential 
source of errors, which arises from the system-bath interactions, can easily be overlooked 
because it is eliminated as soon as semi-classical conditions are imposed. 



IV. DISCUSSION 

The results derived above provide a simple simulation algorithm for the dynamics de- 
scribed by the QCLE. Most often it is the average value of an operator (or correlation 
function) that is of interest. The average value of a quantum operator Bw{X,t) is given by 

B{t) = J dXTi {Bw{X)pw{X,t)) (44) 
= j dX B^\X,t)pi\X), 
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where the trace is taken in the quantum subsystem space. Using Eq. (36) for the time 
evolution of B^'{X,t), the average value may be computed by sampling over the coherent 
state variables and initial density matrix element p^'^{X). 

Our solution for B^' (X, t) has a number of elements in common with other approaches 
that have been devised to simulate nonadiabatic dynamics and it is instructive to make 
comparisons with methods that have been constructed in a similar spirit. 



A. Comparison with partially linearized path integral methods 

First, we draw comparisons between two mixed quantum-classical formalisms: the QCLE 
and partially linearized path integral methods. The formal equivalence between the two 
formalisms was established in a general settin^^ when the subsystem DOF are expressed as 
quantum operators. Therefore, the close resemblance between our solution and that of Huo 
and Cokei® is expected, since they are approximate solutions to the QCLE and a particular 
form of the partially linearized path integral, respectively. However, in view of the derivation 
of our solution presented above, the result in Ref. [20] is not an exact solution of QCLE. 



In our formalism, Hd defined in Eq. (21) contains Vo(-R) = Vb{R) — Tr^/i instead of simply 
the bath potential Vb{R). Recall that the trace term arose from the commutation relation 
for the annihilation and creation operators and the need to use an anti-normal order for 
the product of these operators to evaluate the short-time propagator. If this trace term is 
absent one can show that the solution does not satisfy the differential form of the QCLE. 

The system Hamiltonian, Hiy{X) = Hb{X) + h{R), can be written in an equivalent form 
Hw{X) = Hb{X) + (Trsh{R))/N + h{R), where h{R) is traceless. Since this is an identity, 
the QCLE is independent of the choice of the form which is used in this equation. Our 
solution is also independent of the way the Hamiltonian is written, although the equations 
of motion take a somewhat different form. If the Hamiltonian with the trace removed is 



used in the derivation, the evolution equations have the same structure as is in Eq. (39) 



but Hd{X,z) in Eq. (21) is replaced by HaiX,z) = HoiX) + h^^' zlzy with Hq{X) 



Hb{X) + {Trsh{R))/NJ^When the calculation given in Sec. HI C is repeated with this form 
of the Hamiltonian the QCLE is again obtained, confirming that the different but equivalent 
forms of the Hamiltonian yield the same evolution. 

However, this is not the case when other approximate theories are considered. In partic- 
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ular, it was showiP^ that the choice of Hamiltonian form is crucial in the Poisson Bracket 
Mapping Equation (PBME) approximation to the QCLE (discussed below). When the 
traceless form is used, dynamical instabilities that arise in the course of the evolution can 
be tamed, while if the original form of the Hamiltonian is used the instabilities can lead to 
difficulties. 

The form of the Hamiltonian also affects the nature of the dynamics in the semi-classical 
approach used in Ref. [20]. While the evolution equations in this approach differ from those 



in Eq. (39), the equivalence is restored between the two solutions if the traceless form of 
the Hamiltonian is used. The reason that the partially linearized path integral solution 
depends sensitively on the form of the Hamiltonian is due to the semi-classical approach 
used to solve the dynamics. According to the semi-classical calculation, the dynamics of the 
bath momenta are governed by the force, — | (^dHci{X,z)/dR + dHci{X,z')/dIij, where 
Hci{X,z) = Hh{X) + hci{R,z). The Hamiltonian Hci{X,z) misses the term —Trsh{R) in 
Hci{X, z) in Eq. ( [2l| ) in the current formulation. This extra term is required to restore the 
equivalence between the solution using the original Hamiltonian and that using the traceless 
form of the Hamiltonian. 



B. Comparison with Poisson bracket mapping equation 



Next, we compare the current solution to the PBME approximation to the quantum- 
classical Liouville equationpHSU^ which is obtained from the mapping form of the QCLE by 
dropping an excess coupling termpH. In the case of an isolated subsystem, one can perform 
a change of variables z = {z + z')/2 and A2; = z — z' and show that both the mean, z, and 
the difference, Az, variables follow exactly the same Hamiltonian dynamics, as described in 



Eq. (27) with no R dependence. This implies that if z(0) = A2;(0) then z = Az{t) for all t. 



Since the computation of the time evolution of an operator in the subsystem basis requires 



integration over the entire coherent state phase space, as prescribed in Eq. (38), Az becomes 



a redundant variable. A direct comparison between the two methods can be made if one 
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either integrates out or replaces the integral of A2; by integral of z as follows, 

+ |<^a,a'5m,a^') ■ (45) 

The above identity can be easily proved in a basis that diagonalizes the Hamiltonian, followed 
by transformation of the resulting identity back to the original basis, in the same spirit as 
the computation of the exact coherent state dynamics in Appendix A. 



After properly removing A2;, one can show that Eq. (38) reduces to 



-^^^(t)^;,(t)5,,v + ^5a,a'5^,.') 40(x)5{^^'(X,t), 



= Y.\ dxgx,x'{x)B^\x,t)c,,,{x{t)), (46) 

where x = {q = \/2h^z,p = \/2h^z) and the function^ g\\'{x) = {\mx) {rriyDw and 
c^^>{x) = {d\axf)w represent the Wigner transformation of the outer product of states 
and a pair of annihilation and creation operators, respectively. The last expression in 



Eq. (46) is exactly the evolution of By^'(X,t) in the PBME method. Furthermore, the 
Wigner transformation variables, x, in the PBME method follow the same Hamiltonian dy- 
namics derived above for the mean coordinates of the coherent state variables. Despite the 
very different starting points of the two solutions, this comparison reveals the close relation 
between the dynamics of Wigner transformed coordinates and the mean coordinates in the 
coherent state phase space. Although, this close relation can only be made obvious after 
the effects of difference variables are properly taken into account of and removed (either 



explicitly integrated out or replaced using Eq. (45)). Essentially, the dynamical information 
encoded in the coherent states variables of 2N harmonic oscillators can be merged and be 
encoded in N Wigner transformed coordinates. 

We next comment on the comparison to the PBME method in the presence of a bath. One 
may linearize the bath potential Ve{X, z, z') {Vci{R, z) + ^Az + Vci{R, z) - ^^z)/2 = 
Vci{R, z) such that the dynamics of the bath variables only depends on the mean coordinates 
z. The bath potential linearization allows one to properly remove the difference variables 
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and encode the approximate dynamics in harmonic oscillators. Repeating the same 
calculations and using the coherent-state orthogonality approximation, one can show that 
Eq. (46) still holds in the general mixed quantum-classical setting. 



C. Comparison with semi-classical schemes 

Finally, we compare our results to some semi-classical schemes. The mapping repre- 
sentation consolidates the way subsystem dynamics is handled in the QCLE and in some 
semi-classical schemes. For instance, in the case of a linearized bath potential, the Hamil- 
tonian dynamics prescribed by our solution is also identical to that in the semi-classical 
path integral approach of Stock and Thosd^^'^ as well as the linearized semi-classical initial 
value representation (LSC-IVR) of Millei'^'^^'^. Furthermore, the full version of the current 
solution also handles the subsystem dynamics in ways similar to the forward-backward semi- 
classical initial value representation (FB-IVR) approache^^^'^^'^ that uses the Herman-Kluk 
propagator. One difference is that the forward and backward trajectories are not linked in 
the present solution. 

Finally, we observe that the classical-like system-bath dynamics prescribed in our solution 
could be similar to that of a mixed semi-classical schem^^ in which the bath DOF and sub- 
system DOF are treated with LSC-IVR and the FB-IVR, respectively. Further investigations 
into the subtle connections between the our solution of the QCLE and other semi-classical 
schemes might inspire further developments in nonadiabatic quantum dynamics. 
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APPENDIX A: EXACT EVOLUTION OF COHERENT STATES 

We restrict this analysis to real- valued, symmetric, quadratic Hamiltonian operators, hm, 
which are the only type of Hamiltonian encountered in the mapping formalism. It is always 
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possible to diagonalize such a Hamiltonian matrix, to obtain 

A,A' A, A' ti 

At \ A / V A' 

= Y.KKK^ht, (47) 

where the operators o^^ and are defined in the second fine of the equation. We use the 
superscript d to emphasize that the Hamiltonian is now put in the diagonal form with respect 
to operators and foj^. Since the Hamiltonian is real and symmetric, the matrix M is an 
orthogonal matrix. 

With respect to the new operators and we define the coherent state \y) by 

K\y)=y,\y), {y\hl = {y\y;, (48) 

where 2/^= -^{q^ + ip^). 

Consider time evolution of the coherent state \y) with degrees of freedom, 

e-^ \y) = ®;Li <^ e-l^^l' V ^ |m)^ \ , 



L m=0 J 



= bW) , (49) 

where yv{t) = 2/i/(0)e~^*. In this calculation we used the expansion of a coherent state in 
terms of a complete set of harmonic oscillator states: 



oo 

V" 



7/ 

) = (50) 



Equation (49) implies the equation of motion. 



where /ij?; = -'■^ substitute in the variables g and p into the equation of motion 

for y then we get the usual Hamilton's equation for q and p. 
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Next, we prove that Hamilton's equation is invariant under the hnear transformation, 
Vt^ = Ea -^Ja^a, and y* = J2x ^a-^a^- This proceeds as follows: 



// /I, A' 

A' A 



APPENDIX B: MATRIX ELEMENTS OF THE UNITARY EVOLUTION 
OPERATOR IN THE SINGLE EXCITATION SUBSPACE 

In this Appendix we evaluate matrix elements of the form {m\ \ e"^'*'" \m\i), where hm is 
still the real- valued and symmetric Hamiltonian considered in Appendix A. We evaluate this 
matrix element in two ways: directly and also in terms of matrix elements {fnx\ e^s^™* \frL\i) 
via a linear transformation. The state \m\) = |0i . . . Ij^ . . . Ojv) is an N-harmonic-oscillator 
state with a single excitation on the A-th oscillator, hjb^bx. 

First, we prove that \m.x) — X^^Mj;^ 1"^/*)- This is straightforward since \mx) — a\ |0) 
and \m^) — 6ji |0), where |0) is the common ground state. Therefore, the two states are 
related by the orthogonal matrix M, which was used to establish the hnear transformation 
between a\ and The evaluation proceeds as follows: 



I 



(^^^AV'Mj, {m,,\e-rh'^ \y) {y\ m,) 
dec 

^^-^MavMJ, {m,,\y{t)) {y\ fn,) (53) 



j^Mx,Mi^,x'yAt)y>-^ 

dx , , _| 12 [ dx , , 
(^-A'W^Ae I I = y J^Mt)zxe 



where dx — dqdp and dx = dqdp. To obtain the above result we used the relation Zfj, — 
^vnVf — ^nvyv to re-express the y variables in terms of z variables and employed the volume 
element transformation, dx — dx\<\ei[dya/ dz^]\ — |detM| = dx, since |detM| = 1. 
Since the y{t) variables satisfy Hamilton's equations, \y{t)\'^ — \y\'^- Finally, we note that 
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\y? = T^xVlV^ = T.x,f,y M^'xMxi^zl,z^, = Em^m^m = completing the results needed to 



obtain Eq. (53) 



Next, we compute e n^'^\z) directly, where \z) is defined by dx \z) = zx \z). To carry out 



this calculation we reconsider Eq. (53) 



{mx'\e '^^-^^\mx) 
dx 



dx 



N 



{2nh) 



N 



[mx'\e ft''™* 



{mx'\e ft^™* \z) {z \mx) 
z) 2;^e"^|2|^ 



(54) 



Comparing the last lines of Eqs. (53) and (54), we see that (mv|e r"™* 1^;) = Zx'{t)e 2 



\\<t)? 



(mv| z{t)). Since the identities hold for all possible {mx'\ and 1^;), we can identify e \z) 
\z{t)). 



APPENDIX C: EFFECTIVE LIOUVILLE OPERATOR 



Below, we prove the identity in Eq. (29) that relates the forward and backward bath 



propagators to the effective Liouville operator: 



j=Q k=0 V"' V / 



j=0 ^' k=0 {p} 



f <^(i(£ + 0)' iw-(A-) 

i=o -J- ^ / 



= e^^=(^'^'^')Miy (55) 

In these expressions we used the shorthand notations, C=C {X,x) and C'=C (X, x'), and 
the definition of S in going from the third to fourth lines. The sum on {p} denotes a sum 
over all permutations. 
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APPENDIX D: DIFFERENTIAL FORM 



The term zeroth order in r in Eq. (42) is easily computed by performing the integrals 
over z{t) and z'{t) and the result is simply {X,t). The first term of order r, which we 
call Ji, involves time derivatives coherent state variables. Using the equations of motion for 
the coherent state variables and performing the integrals we find 



h='-{\\[h,Bw{X,t)]\X') 

= '-{X\[Hw,Bw{X,t)]\X'), (56) 



which is the first term in the QCL operator in Eq. ([2]). 

Next, we consider the first-order term involving the evolution of the spatial coordinates 



of the bath as given by the effective Liouville operator. Inserting its definition in Eq. (30), 
iCe{X, z, z') = ■ -§ji — ^^"^qr'^ ■ ^) the evaluation of the term involving ■ 
is straightforward since it does not contain the coherent state variables. Performing the 
integrals over these variables yields I2 = ■ The remaining terms require more 

attention since they involves the force acting on the bath variables, which depends on the 
effective potential where Ve{X, z, z') = {Vci{R, z) + Vci{R, z'))/2. Denoting this contribution 
Is, we have 



~Hj d^4t)d^z'{t)(f){z)(j){z') 
-J2f dh{t)dh'{t)(f){z)(j){z') 



xz,it)z'*,it)z*(t)z'it) 



dP 



dVojR) l dV,i{R,z) l dV4R,z') 
dR 2 dR 2 dR 



dB^ dVojR) 
dP OR 



+ h, + I: 



32' 



(57) 
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The integral may be evaluated as follows: 

^3i = -^E / dh{t)dh\t)cl>{z)cl>{z') 

x^>^W^x'W^,,WZh'W-qP — 

= -^E/ d'zmz)zx{t)z;{t) 



ixaa 



dP dR 



dP dR 

ixaa' ' 

[\zxit)\^\z^{t)\^5aa'S^.x{l - S^X) 
+ \zx{t)\^Saa'SaxSlxx] 

Performing the z integrals we find 



'3i 



olZ^ PlR Pip 



OR dP ^ dR dP 

dR dP ) 

- "2 ^'^i + ^^^^ 

1 fdhdBw , aTrs/iafivF\ I w\ /rox 

- "2 la^^ + ^ ■ ^^^^ 

In writing the last line of this equation we used the fact that the subsystem Hamiltonian is 
independent of so and be replaced by h. 

Similarly, the Is^ integral can be evaluated to give, 

1 (dB,vdh dBwdTTsh] 

Recall that Vo{R) = -Tr,/i so that ^^gf ^ = - Given these 

results the entire Is integral is 

h = -I ((A| {Hw,Bw} {Bw,Hw} |A')) , (60) 
where the Trg/i terms arising from the Vq, Is^ and I^^ canceled. 
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FIG. 1. (Color online) Schematic diagram of the time evolution of the bath coordinates X = {R, P) 
(the green line), the forward coherent state coordinates x = {q,p) (the blue line), and the backward 
coherent state coordinates x' = {q',p') (the red line). The vertical axis denotes the time. At 
each time step i, the classical Hamiltonians Hci{Xi, Xi) and Hci{Xi, x[) are parametrized with the 
updated coordinates. The wiggly, orange lines represent the direct coupling between the evolutions 
of different sets of phase space coordinates under the influence of the classical Hamiltonians. As 
shown, the two sets of coherent state variables are only coupled via the bath coordinates. 
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